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Multifractal Scaling of Thermally- Activated Rupture Processes 
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We propose a "multifractal stress activation" model combining thermally activated rupture and 
long memory stress relaxation, which predicts that seismic decay rates after mainshocks follow the 
Omori law ~ 1 /t'' with exponents p linearly increasing with the magnitude Ml of the mainshock and 
the inverse temperature. We carefully test this prediction on earthquake sequences in the Southern 
California Earthquake catalog: we find power law relaxations of seismic sequences triggered by 
mainshocks with exponents p increasing with the mainshock magnitude by approximately 0.1 — 0.15 
for each magnitude unit increase, from p{Ml = 3) ~ 0.6 to p{Ml — 7) ~ 1.1, in good agreement 
with the prediction of the multifractal model. 
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Parisi and Frisch 1] and Halsey et al. Q have intro- 
duced the extended concept of scale invariance, called 
multifractality, motivated by hydrodynamic turbulence 
and fractal growth aggregates respectively. Use of the 
multifractal spectrum as a metric to characterize com- 
plex systems is now routinely used in many fields, in- 
cluding seismology to describe the hierarchical structure 
in space and time of earthquakes and faults (see for in- 
stance 1^ Q 111). However, the origin of multifractal- 
ity is rarely identified. This is certainly true for earth- 
quakes for which the possible existence of multifractal- 
ity is under scrutiny due to limited and corrupted data 
sets leading to biases and its origin a matter of de- 
bate: fractal growth processes 0, self-organized critical- 
ity 3 or hierarchical cascades of stresses are among 
the physical scenarios proposed to lead to multifractal- 
ity in fault and earthquake patterns. Here, we propose a 
physically-based "multifractal stress activation" model of 
earthquake interaction and triggering based on two sim- 
ple ingredients: (i) a seismic rupture results from ther- 
mally activated processes giving an exponential depen- 
dence on the local stress; (ii) the stress relaxation has a 
long memory. The interplay between these two physical 
processes are shown to lead to a multifractal organization 
of seismicity, which we observe quantitatively in real cat- 
alogs. 

Thermal activation is relevant in all previously pro- 
posed physical processes underlying earthquakes: creep 
rupture, stress corrosion and state-and-velocity depen- 
dent friction. We model seismic activity A(r, t) at po- 
sition r and time t as the occurence of frictional slid- 
ing events and/or fault ruptures that are thermally ac- 
tivated processes facilitated by the applied stress field: 
A(r, t) ~ exp [—PE{f, t)], where /3 is the inverse temper- 
ature and the energy barrier £^(r, t) for rupture can be 
written as the sum of a contribution EQ{f) characterizing 
the material and of a term linearly decreasing with the 



locally applied stress S(f, t): E{f, t) = Ea{r) - Vi:{r, t). 
y is a constant which has the dimension of a volume and 
S(r, t) is the total stress at position r and time t. The 
decrease of the energy barrier E{r, t) as a function of the 
applied stress £(r, t) embodies the various physical pro- 
cesses aiding rupture activation under stress. In addition, 
there are many evidences for a stress-controlled earth- 
quake activation process, suggesting that earthquakes 
trigger earthquakes directly and indirectly via dynami- 
cal and static stress transfers. Visco-elastic models of 
stress relaxation can account for the short-term relax- 
ation processes of the strain measured by geodetic meth- 
ods but, over long time scales, it is necessary to take 
into account the presence and geometry of lower crustal 
and mantle shear zones, which lead to slower decay- 
ing relaxation rates. We thus write the stress S(r, t) 
at position r and time t as the sum of contributions 
from all past events at earlier times t < t and posi- 
tions r ': S(f,i) = Sfar field (r, t) + j'^^ J dN[dr ' x 
dT]Aa{r ' ,T)g(r — f' ,t — T). A given past event at (r ',t) 
contributes to the stress at (r, t) by its stress drop am- 
plitude A(T(r ',r) which is transfered in space and time 
via the stress kernel (or Green function) g{r — ,t — r), 
taking into account both time relaxation and spatial ge- 
ometrical decay. The term dN[dr ' x dr] is the number 
of events in the volume df ' that occurred between t and 
T -f dr. 

In this letter, we restrict our analysis to the time do- 
main. For this, we assume for simplicity that g{r,t) is 
separable as g{r,t) — f{r) x h{t). This obtains 



\{f,t) = Xtcc{r,t) exp 



/? / dr s{r, T)h{t ~ t) 



(1) 



where s(r, t) = J df ' Acr(r ',r) fif—f') is the effective 
source at time r at point r resulting from all events occur- 
ring in the spatial domain at the same time r. Atcc(^, i) 
is the spontaneous seismicity rate in absence of stress 



2 



triggering by other earthquakes and accounts for the tec- 
tonic foading (far field stress), which may in general be 
non-homogeneous in space and perhaps depends on time. 
Since expression JQ) is defined for any r, we drop the ref- 
erence to r without loss of generality. 

To go further, we specify the distribution P{s) of 
stress sources and the memory kernel h{t). On the ba- 
sis of theoretical calculations, simulations and measure- 
ments of rotations of earthquake focal mechanisms, Ka- 
gan has suggested that P{s) should follow a symmetric 
Cauchy distribution. To capture in a phenomenological 
way the extended nature and complexity of earthquake 
ruptures, we use a more general power law distribution 
P{s) ^ C/|As|^+^, which generalizes the Cauchy case 
/i = 1. To account for the slower-than-exponential stress 
relaxation processes discussed above, we postulate that 
^i't) = (t+l|i+tf for < f < T, which is of the Omori 
form with the usual small time-scale cut-off c. To ensure 
convergence of the correlation function of deterministic 
processes with memory governed by h{t) for any possible 
values of 9, we truncate the power law at some large time 
T, which we call the "integral time scale:" it is the largest 
time scale up to which the memory of a past event sur- 
vives. T can thus be interpreted as the effective Maxwell 
time of the relaxation process. The time dependence of 
h{t) is an effective description of the relaxation of stress 
due to microscopic processes such as dislocation motion, 
stress corrosion and hydrolytic weakening which obeys 
an Omori-like power law. 

In summary, our model reads (in discretized form) 



get E[e' 



.uj(t) 



A(i) = Atoc e 



with the stress sources s(ti) distributed according to a 
power law P{s) with exponent fi and h{t) having a power 
law memory. 

We now derive our novel prediction for Omori's law 
quantifying the decay of seismic activity after a "main- 
shock" occurring at the origin of time. This amounts to 
determining the typical time dependence of the seismic 
rate X{t) conditioned on a value Am realized at t = 
which is larger than average. This formulation is due to 
the fact that a mainshock of magnitude M induces a local 
burst of seismic activity proportional to K 10"^, where 
K and a are two positive constants . Since the stress 
sources are non-Gaussian but power law distributed, 
their average and variance may not be defined. Rather 
than calculating the conditional expectation of A(i), a 
typical measure of conditional seismicity rate can be de- 
fined at any quantile level q by the probability Pr[A(t) > 
A^IAm] that the rate A(i) be larger than the quantile 
Xq conditioned on the fact that the seismic rate was at 
some given value Xm at time 0: Pr[A(t) > AgjAjv/] = 
p^[g/3^(t) > ^\u;m] = PrHt) > (l//3)ln(A_) |^^,]. 
For Gaussian sources, lu is normally distributed and we 



Using lO, this 



ujm] = exp PE[uj{t)\ujMi 1 2 

where E[^(t)|^M] - %t[tT^ ^ 
would provide a closed formed expression for the Omori 
law describing the relaxation of the conditional rate 
E[A(t)|AM]- The physical meaning of this result is that 
one can write a linear regression u;(t) = ^(t)ujM+£, where 
7(t) is a non-random factor and e is a centered Gaus- 
sian noise with zero correlation with ujm ■ This equation 
writes that the best predictor of uj given lum is jlum: i-e., 
F,[uj{t)\ujM] = JUJM with 7 = ^°vL" [ !.ay"'^ • power law 
stress sources, we use the insight that the natural gener- 
alization of the variance for power laws p{x) « C/x^'^'^ 
with infinite variance (i.e., with /i < 2) is the scale pa- 
rameter C (see Chap. 4 of ^3|)- In the power law case, 
due to the linear form of to in we can still write 
Lu(t) = ^{t)ujM + e but with uj{t),ujM and e being power 
law distributed random variables with the same expo- 
nent /i and with scale factors equal respectively to C^^ 
(for UJ and um) and C^. The key idea is that 7 can 
be determined by forming the random variable defined 
as the product ujlum = + ^^m- It is straightfor- 

ward to show that the distribution of u)U>m consists of 
two main contributions, (i) a dominant power law with 



exponent and scale factor C^. 



= 7^/^ C^, and 



(ii) a sub-dominant power law with exponent /i (with 
a logarithmic correction) and scale factor C^jC^. This 
has the following practical implication: if one measures 
or calculates the leading power law decay of lu x lum, 
the measure of its scale factor gives access to the pa- 
rameter 7 through the expression 7(f) = {C^^i^,^ /C^)i^ . 
where the time dependence of ^(t) comes from that of 



For /i 



2, we recover the Gaussian result 



with the correspondence G^^ = Var[a;] and C^jujm = 
Cov[[j(t), WAf]. Using we then form the product 

u:{t)ujM = E,; I u<tY.j I t,<o ''^fe) Ht-U) H-tj), 

where the s's are random variables with power law tail 
with exponent /i. Then, using standard calculations (see 
Chap. 4 of 11]), the terms in the double sum that con- 
tribute to the leading asymptotic power law tail with ex- 
ponent /i/2 correspond to the diagonal terms i = j, while 
all the other terms contribute to the sub-leading power 
law tail with exponent /i with logarithmic corrections. 
This gives the expression of the scale factor ci,(i{^'' of 

the dominating power law with exponent fj,/2 and finally 

2 

yields 7 = (^J2i \ uKoiKt - ti)h{-ti)]^y\ in discrete 
form and 



7(t) 



1 



dy 



1 



1 



c/t 



(y-l-l)™ 



(3) 

in continuous form where m = (1 -I- 9)fj,/2. The discrete 
time step At converting the discrete into the continu- 
ous sum is the average time interval between two events 
before a mainshock. 
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We thus obtain Pr[w(t) > y\ujM] = Pr[7t^A/ + e > 
y\ujM] = Pr[e > y - "tu}M\^M] = F{y - j{t)LL>M), 
where F(e) is the complementary cumulative distribu- 
tion of e. Putting these results in this leads to 

Pr[A(t) > A, I Am] = In (^) -jit)u;M). The 

typical time evolution of the seismicity rate X{t) condi- 
tioned on the rate Am at time is thus given by fixing the 
quantile probability to some level Pr[A(i) > Aq|Aj\/] = q, 
leading to 



Xg{t) = Ag At 



(4) 



where Aq = exp (/3f'~^(q)) . The time-dependence of the 
seismic decay rate requires the determination of the time- 
dependence of j{t) given by ©). We now show that, for 
a rather broad range of values of the exponents fi and 9 
defining the model, Xq{t) is approximately given by 
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p{M) = al3M + 6/3 , 



(5) 



where a > and M is the mainshock magnitude. 

Consider first the case 2m — ^{1 + 6) = 1, such that 
the exponent to = (1 -I- 0)^/2 defined in equal to 



1/2. Then, 



m/2 



dt 
dy" 



is close to 



(T+c) 



1/2 



(T + C-t)l/2 (t + c)l/2 

l/t, and thus 7^/2 (t) 



showing that 

constanti — constant2 \n{t/T) which, for not too small nor 
too large <'s and for constanti < constant2, gives « 
constant'^ — constantj x \n{t/T). This yields ((SJ). Typi- 
cally, the power law behavior is observed over more than 
two decades in time, which is comparable to empirical 
observations, as verified by direct numerical integration 
of (O. Then, expression Q leads to (|3J| using the fact 
that LUM oc In(AM) oc Iti{K 10"*^ = a In 10 M -h InX, 
i.e., ujm is linearly related to the magnitude M. The fact 
that ^{t) is asymptotically exactly logarithmic in time 
for 2to = /x(l + 9) ^ 1 and thus that the seismic rate 
X{t) is an Omori power law can be recovered from a dif- 
ferent construction motivated by multiplicative cascades 
introduced in turbulence ^3 ■ This case covers the exact 
multifractal random walk model [l^ , which corresponds 
asymtotically to 9 ~ —1/2 and fj, = 2. This continu- 
ous dependence of the exponent p{M) has actually been 
documented empirically in this case in another context 
of aftershock decay following shocks in financial markets 
[1^ . For 2m = /i(l -1-6*) 7^ 1, one can often observe an 
approximate linear decay of "f(t) as a function of In i, over 
two to three order of magnitudes in time in the decaying 
part, all the more so, the closer to is to 1/2, also leading 
to ©. 

We now show that this prediction is verified in 
the Southern Californian earthquakes catalog with re- 
vised magnitudes (available from the Southern California 
Earthquake Center). The details of our analysis is given 
elsewhere jTB ] and we summarize the main results. In 
order to improve the statistical significance and to test 



for the stability of our analysis, we analyzed four differ- 
ent sub-catalogs: 1932 — 2003 for magnitude Ml > 3 
(17,934 events), 1975 - 2003 for Ml > 2.5 (36,614 
events), 1992 - 2003 for Ml > 2 (54,990 events), and 
1994 - 2003 for Ml > 1.5 (86, 228 events). We consider 
all events in a given sub-catalog and discriminate between 
mainshocks and triggered events ("aftershocks"). Main- 
shocks are determined by using two different decluster- 
ing methods described below. Once the mainshocks are 
determined, triggered events are defined as those events 
following a mainshock, which belong to a certain space- 
time neighborhood of it. In order to test for the predicted 
dependence of the p-value as a function of magnitude, 
we bin the mainshock magnitudes in intervals [1.5; 2], 
[2; 2.5], [2.5; 3], and so on up to [7; 7.5]. In each main- 
shock magnitude interval [Mi ; M2] , we consider all trig- 
gered sequences emanating from mainshocks with mag- 
nitude in this interval and stacked them to a common 
origin of time. The resulting function is fitted using the 
modified Omori law N{t) = B + pq^^, where B is a 
positive parameter introduced to account for the back- 
ground seismicity assumed to be superimposed over the 
genuine triggered sequences. The time shift c ensures the 
regularization of the seismic rate at t = 0. 

The first declustering method is essentially the same 
as defined in ^3] ■ every event in the catalog is defined as 
a mainshock if it has not been preceded by an event with 
larger magnitude within a fixed space-time window Txd, 
with T = 1 year and d = 50 km. Looking for events trig- 
gered by this mainshock, we define another space-time 
window following it. The time dimension of the window 
is also set to 1 year, whereas the space dimension depends 
on the rupture length of the main event. This spatial win- 
dow is chosen as a circle of radius equal to the mainshock 
rupture length L = io~2-57+o.6Mi, ^ which is the average 
relationship between L and magnitude Ml for California 
|l6l] . If any event falls within this space-time window, it 
is considered as triggered by the main event. We have 
also checked the stability of the results by considering a 
spatial neighborhood of radius 2L rather than L for the 
triggered events. The second declustering method is the 
same as the first one, except for one element: the space 
window used for qualifying a mainshock is not fixed to 
d = 50km but is chosen to adapt to the size of the rupture 
lengths L{M,) given by L = iq-2.57+o.6Ml of events 
of all possible magnitudes ML{i) preceding this potential 
mainshock. 

Figure ^ shows sets of typical seismic decay rates of 
stacked sequences for several magnitude intervals of the 
mainshocks, for the period from 1932 to 2003 when using 
the first declustering technique, with mainshock magni- 
tudes above Ml = 1.5. Very similar plots are obtained 
for different time periods, with the second declustering 
method and by varying the size from L to 2L of the spa- 
tial domain over which the triggered sequences are se- 
lected For large mainshock magnitudes, the roll-off 
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FIG. 1: Normalized seismic decay rates of stacked sequences 
for several magnitude intervals of the mainshocks, for the pe- 
riod from 1932 to 2003 when using the second dedustering 
technique. 



at small times is due to the observational saturation and 
short-time lack of completeness of triggered sequences. 
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ETAS model provides a particular strong null hypothesis 
as it rationalizes most of the phenomenological statisti- 
cal properties of earthquake catalogs 17]. By construc- 
tion, synthetic catalogs generated with the ETAS model 
should exhibit Omori laws with magnitude-independent 
exponents. Applying our procedure to such synthetic 
catalogs allows us to investigate whether the magnitude- 
dependence of the p-value reported above could result 
from some bias introduced by our analysis rather than 
being a genuine property of earthquake catalogs. We 
verify that p{M) obtained by our procedure is a con- 
stant independent of M equal to the input value used in 
the generation of the synthetic catalog (l^ . 

Let us conclude by offering an intuitive explanation 
of © using the properties of multifractal spectra. The 
temporal evolution of seismicity in a fixed spatial domain 
defines a statistically stationary measure on the tempo- 
ral axis, the measure determining the rate of earthquakes 
at any possible instant. An Omori sequence with expo- 
nent p corresponds to a singularity (to the right) equal to 
1 — p (logarithmic for p = 1). A large earthquake triggers 
a strong burst of seismicity, giving rise to a strong singu- 
larity. For the relation a ~ 1 — p to be consistent with 
the multifractal description, a large earthquake must be 
associated with a strong singularity, a small a, hence a 
large p. Reciprocally, small moment orders q select weak 
seismic sequences, which are thus associated with small 
local mainshocks. Small q's are associated with large a's, 
hence small p's. By a similar argument in the space do- 
main, the exponent of the spatial decay of the seismic 
rate induced by a mainshock of magnitude should 
increase with M^. Thus, in this view, the ETAS model 
is nothing but the mono-fractal approximation of a more 
general multifractal description of seismicity. 

This work was partially supported by NSF-EAR02- 
30429 and by the Southern California Earthquake Center 
(SCEC). 



FIG. 2: p- values of the Omori law ~ obtained by the 
procedure described in the text for mainshocks (defined using 
the second dedustering algorithm) as a function of the main 
events' magnitude, for the difi'erent sub-catalogs of lifespans 
given in the inset. 

Figure [3 shows the fitted p- values as a function of the 
magnitude of the mainshocks for each of the four sub- 
catalogs. We use a standard least-square fit of the seismic 
rate as a function of time with a weight proportional to t 
for each bin to balance their relative importance. We also 
take into account the possible presence of a background 
term. We have also performed maximum likelihood esti- 
mations of the exponent p, confirming the results shown 
in Fig. 12 . To test the reliability and robustness of 
our results, we have simulated synthetic catalogs with the 
ETAS model with known statistical properties following 
exactly the same procedure as for the real catalogs. The 



* Electronic address: Isornet t e@moho .ess ..ucla.edu | 
[1] Parisi, G. and Frisch, U., in Proc. Int. School Enrico 
Fermi, eds. M. Ghil et al. (North Holland, Amsterdam, 
1985). 

[2] Halsey, T.C. et al., Phys. Rev. A 33, 1141 (1986). 
[3] Godano C. et al., Geophys. J. Int. 136, 99 (1999). 
[4] Main, I., Rev. Geophys. 34, 433 (1996). 
[5] Ouillon G. et al., J. Geophyss Res. 101, 5477 (1996). 
[6] Ouillon, G. and Sornette D., Geophys. Res. Lett. 23, 3409 
(1996). 

[7] Sornette, A. et al., J. Geophys. res. 98, 12111 (1993). 
[8] Rodkin, M.V., Izvestiya-Physics Sol. Earth. 37, 663 
(2001). 

[9] Kagan, Y.Y., Nonlinear Proc. Geophys. 1, 171 (1994). 
[10] Helmstetter, A., Phys. Rev. Lett. 91, 058501 (2003). 
[11] Sornette, D., Critical Phenomena in Natural Sciences, 

2nd ed. (Springer, Heidelberg, 2004). 
[12] Schmitt, F. and Marsan, D., Eur. Phys. J. B 20, 3 (2001). 
[13] Muzy, J.-F. and E. Bacri, Phys. Rev. E 66, 056121 (2002). 
[14] Sornette, D. et al., Risk 16, 67 (2003). 



5 



[15] Ouillon, G. and D. Sornette, cond-mat/0407208 [17] Helmstetter, A. and D. Sornette, J. Geophys. Res. 108, 

[16] Wells, D.L. and Coppersmith, K.J., Bull. Seism. Soc. 10.1029, 2003. 

Am. 84, 974 (1994). 



